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Abstract 

In this paper, we develop a multiscale finite element method for solving flows in fractured media. Our 
approach is based on Generalized Multiscale Finite Element Method (GMsFEM), where we represent 
the fracture effects on a coarse grid via multiscale basis functions. These multiscale basis functions are 
constructed in the offline stage via local spectral problems following GMsEEM. To represent the fractures 
on the fine grid, we consider two approaches (1) Discrete Fracture Model (DEM) (2) Embedded Fracture 
Model (EFM) and their combination. In DEM, the fractures are resolved via the fine grid, while in EFM 
the fracture and the fine grid block interaction is represented as a source term. In the proposed multiscale 
method, additional multiscale basis functions are used to represent the long fractures, while short-size 
fractures are collectively represented by a single basis functions. The procedure is automatically done via 
local spectral problems. In this regard, our approach shares common concepts with several approaches 
proposed in the literature as we discuss. Numerical results are presented where we demonstrate how one 
can adaptively add basis functions in the regions of interest based on error indicators. We also discuss 
the use of randomized snapshots ([2]) which reduces the offline computational cost. 


1 Introduction 

Many porous media flow and transport processes are dominated by the presence of fractures. The fractures 
present high conductivity conduits and their effects need to be captured accurately in the simulations. 
Because the fractures have very small width, they are often represented as zero thickness objects in numerical 
simulations. Because there are many fractures of different lengths (connected or not connected), the problem 
is inherently multiscale and efficient multiscale techniques are needed to represent their effects. 

In this paper, we present an application of Generalized Multiscale Finite Element Method for flows in 
fractured media. The study of flow in fractured media on the fine grid has been conducted in numer¬ 
ous papers. These models applied various fracture models, such as the Discrete Fracture Model (DFM), 
Embedded Fracture Model (EFM), the single-permeability model, and the multiple-permeability models 
([111 [niioiiiTi nans]). Though these approaches are designed for fine-scale simulations, a number of these 
approaches represent the fractures at a macroscopic level. For example, multiple-permeability models repre¬ 
sent the network of connected fractures macroscopically by introducing several permeabilities in each block. 
The EFM ([191 UHl Ull) models the interaction of fractures with the fine-grid blocks separately for each 
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block. These approaches can be generalized by incorporating the interaction of fractures and permeability 
heterogeneities locally which can lead to efficient upscaling techniques. Some general upscaling techniques 
for flows in fractured network are presented in 13 HB, where the authors introduce the fracture and matrix 
interaction parameters hierarchically. 

In recent papers [13], several multiscale approaches are proposed for representing the fracture effects. 
These approaches share common concepts with the methods that we discuss here in a sense that they add 
new degrees of freedom to represent the fractures on a coarse grid. The main difference is that our approaches 
use local spectral problems accompanied by adaptivity to detect the region where to add new basis functions. 
In this regard, the procedure of finding multiscale basis functions and the enrichment procedure is different 
from existing techniques. 

In this paper, we use a general multiscale finite element framework, GMsFEM. GMsFEM is a flexible 
general framework that generalizes the Multiscale Finite Element Method (MsFEM) ([l5]) by systematically 
enriching the coarse spaces and taking into account small scale information and complex input spaces. In this 
work, we use Discrete Fracture Model (DEM) ( [2Tl[T4l [24l[T6]) for the simulation of the fine-scale problem 
as well as in the construction of spectral problem for the GMsFEM. In GMsFEM approach, as in many 
multiscale model reduction techniques, which divides the computation into two stages: the offline stage and 
the online stage. In the offline stage, a small dimensional space is constructed that can be used in the online 
stage to construct multiscale basis functions. These multiscale basis functions can be re-used for any input 
parameter to solve the problem on a coarse grid. The main idea behind the construction of offline and online 
spaces is the selection of local spectral problems and the selection of the snapshot space. In [7], we propose 
several general strategies. The main idea of this paper is to combine GMsFEM with DEM and/or EFM to 
solve the flow problem in fractured media. We present a multiscale basis construction, adaptivity, and the 
use of randomized snapshots to reduce the computational cost. 

Our approaches share some common concepts with hierarchical fracture modeling proposed by Lee et 
al. [m, where the main idea is to homogenize small-length fractures (with the length smaller than the 
coarse block), while represent the large-length fractures. Our methods via local multiscale basis functions 
homogenizes small-size fractures (by lumping their effect into a single-per-node multiscale basis function) 
and represent the long-size fractures in each coarse block. A global coupling, such as finite element in this 
case, recover an accurate representation of long size fractures by coupling all the information together. 

We present several numerical examples to illustrate the performance of the proposed approach. We 
consider the fine-scale fracture representations using DEM, EFM as well as coupled DEM and EFM, where 
the shorter fractures are represented by DEM and the longer ones EFM. All the numerical results indicate 
that the proposed GMsFEM is robust and accurate. Besides, we test the performance of adaptive algorithm 
in m for the problem in this paper. To reduce the computational cost, we propose the use of randomized 
snapshots where only a few snapshots are constructed and used to construct multiscale basis functions. We 
would like to emphasize that our goal is to develop and show the performance of GMsFEM for flows in 
fractured media and we do not make any comparisons between different hue-scale fracture models. 

The rest of the paper is organized as follows. In Section we present preliminaries. The hne-scale 
fracture modeling techniques are briehy reviewed in Section |3.1[ The construction of the coarse spaces for 
the GMsFEM is displayed in Section In Section numerical results for several representative examples 
are showed. Finally, we conclude our paper with some remarks in Section 

2 Preliminaries 

In this paper, we study high-contrast how problem 

— diY (^k{x) u) = f in D, (1) 

where k{x) has fractures with high values and low values in the matrix (see Figurefor illustration). To 
discretize Q, we introduce the notions of hne and coarse grids. Let be a usual conforming partition 
of the computational domain D into hnite elements (triangles, quadrilaterals, tetrahedra, etc.). We denote 
this partition as the coarse grid and assume that each coarse element is partitioned into a connected union 
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of fine-grid blocks. The fine grid partition will be denoted by T^, and is by definition a refinement of the 
coarse grid T^. We use {xi}fLi (where N denotes the number of coarse nodes) to denote the vertices of the 
coarse mesh and define the neighborhood of the node Xi by 

oji = [j{Kj e r^; Xi€Kj}. (2) 

See Figure for an illustration of neighborhoods and elements subordinated to the coarse discretization. We 
emphasize the use of uji to denote a coarse neighborhood, and K to denote a coarse element throughout the 
paper. 
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Figure 1: Illustration of a coarse neighborhood and coarse element 



Figure 2: A heterogeneous fracture field. 

Next, we briefly outline the GMsFEM. We will consider the continuous Galerkin (CG) formulation and 
signify uji as the support of basis functions. We denote the basis functions by which is supported in 
and the index k represents the numbering of these basis functions. In turn, the CG solution will be sought 
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as Ums{x) = Oiice the basis functions are identified (see the next section), the CG global 

coupling is given through the variational form 


a('Ums, v) = (/, v), for all v G Kff, (3) 

where Vqs is used to denote the space spanned by those basis functions and a(-, •), (/, •) are bilinear (or linear) 

form corresponding to (1) defined by a{u^ '^) = / ki{x)Vu • Vu dx, and (/, v)= fv dx. We also note that 

one can use discontinuous Galerkin formulation (see e.g., [6l[3l[9]) to couple multiscale basis functions defined 
on K. 

Let V be the conforming finite element space with respect to the fine-scale partition T^. We assume 
u G V is the fine-scale solution satisfying 


a{u,v) = {f,v), veV. (4) 

In the next section, we will introduce fine-grid discretizations for fractures. As we have mentioned above 
that the aperture of the fracture is very thin, and the fracture permeability is high. 

3 Discretizing fractures on the fine grid 

3.1 Discrete Fracture Model (DFM) 

Our first approach is based on representing the fractures on the fine grid as the edges of finite element 
mesh. This allows meshing the fractures more accurately; however, it can be expensive when the fracture 
distribution is irregular. Following a standard convention, we call “matrix” the region that excludes all the 
fractures. In the DFM (see e.g., [T]) as implemented in this paper, we assume that the permeability does 
not vary along the fractures, which coincide with the edges (or faces) of finite element mesh. The fracture 
aperture is taken into account in the finite element discretization as a lower dimensional lines or surfaces. 
More precisely, the discretization of the fractures on low dimensional surfaces are added to finite element 
discretization of the matrix system. Below, we demonstrate this main idea in our two dimensional example. 

Based on the assumptions above on the DFM, the fractures can be treated as one dimensional. One¬ 
dimensional element is introduced in addition to the two-dimensional element for the discretization of the 
matrix. Thus the Eqn. 0 will be discretized in two-dimensional form for the matrix and in one-dimensional 
form for the fractures. The whole domain D can be represented by 

D = D^U(UiDfrac,i), (5) 

where the subscript m and frac represent the matrix and the fracture regions, respectively and i refers to 
ith fracture. Note that is a two-dimensional domain and Dfrac,i is a one-dimensional domain. We write 
the finite element discretization corresponding to Eqn. Q as (see also Eigurej^for illustration) 

[ k{x)Vu ■Vvdx= [ k{x)Vu -Vvdx + Y^ f k{x)Vu ■Vvdx= [ fv dx. (6) 

D J Dm ^ '^Dfrac,i D 


3.2 Embedded Fracture Model (EFM) 

In this section, we present a fine-discretization of fracture network following EEM proposed by Lee et al. 
(ini)- EEM allows avoiding complex meshing for fracture network. We would like again to emphasize that 
our goal is not to compare DEM with EEM, but rather to develop GMsEEM framework which uses both 
fine-grid discretization techniques. In EEM, long fractures are treated as low dimensional objects and their 
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Figure 3: Illustration of DFM. 


interaction with the matrix is modeled separately. In our two dimensional simulations, we treat the long 
fractures separately as one-dimensional problems. Eqn. § can be written as 

—div(/T:(x) Vi^)-1-/ in (7) 

-diY{K{x) \/u) = / in fractures, (8) 

where and represent the interaction between the long fractures and the rest of the system. Those 
two terms can be calculated following HU- We denote the linear system corresponding to the discrete form 
of the problem above is 



Blf ■ 

mf 




-fra- 



0 


v} 

= 


_ 1 

0 



U^f_ 


jNf_ 


(9) 


Here, Nf denotes the total number of long fractures, is the matrix obtained over the domain excluding 
the long fractures, denotes the interaction of the ith long fracture with the matrix and denotes 

the matrix interaction with the ith long fracture, where i = 1, 2, • • • , Nf. Here, we use EFM to handle long 
fractures, though it can be used for shorter fractures. 


4 GMsFEM 

In this section, we will give a brief description of the GMsEEM for heterogeneous flow problems. More 
details can be found in miH]. In the following, we give a general outline of the GMsEEM. Then, we will 
discuss the multiscale basis construction. 

Offline computations: 

Step 1 Coarse grid generation. 

Step 2 Construction of snapshot space that will be used to compute an offline space. 

Step 3 Construction of a small dimensional offline space by performing dimension reduction in the space of 
local snapshots. 

If the problem has a parameter, one can use the offline space to construct multiscale basis functions in 
the online stage. In the above outline, the offline space can be reused if we change the input of the model. 
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Given the computational domain, a coarse grid can be constructed and local problems are solved on coarse 
neighborhoods to obtain the snapshot spaces. Then, smaller dimensional offline spaces are obtained from 
the snapshot spaces by dimension reduction via some spectral problems. 

4.1 Local basis functions 

We now present the construction of the basis functions and the corresponding spectral problems for obtaining 
a model reduction. 

In the offline computation, we first construct a snapshot space V^nap edich coarse neighborhood uji. 
The snapshot space can be the space of all fine-scale basis functions or the solutions of some local problems 
with various choices of boundary conditions. In [2], we use harmonic extensions to form a snapshot space. 
For fiows in fractured media, the snapshot vectors are computed in the same way except that we use an 
appropriate fine-grid discretization taking into account the fracture distribution. In specific, given a fine-scale 
piecewise linear function 6j{x) defined on duji^ we define by 

— div(/T:(x)V'0^"’®^^^) = 0 in (10) 

where = S^{x) on duji. The local problem is solved taking into consideration the fracture distribution. 

Later, we use randomized boundary conditions to avoid computing the full snapshot vectors and use only a 
few snapshot vectors. In particular, we use DFM for computing all snapshot vectors. For EFM, we follow a 
hierarchical approach and discretize few fractures using the interaction terms. We can also add them into the 
calculations of snapshot vectors; however, in the paper, we consider EFM for handling a few long fractures. 

For brevity of notation we now omit the superscript yet it is assumed throughout this section that the 
offline space computations are localized to respective coarse subdomains. Let k be the number of functions 
in the snapshot space in the region and 

14nap = span{'*/)®“®'P : l<j<li}, 

for each coarse subdomain uJi. Denote 

Rsn., = [rr^,...,rir]. 

In order to construct the offline space we perform a dimension reduction in the local snapshot space 
using an auxiliary spectral decomposition. The analysis in m motivates the following eigenvalue problem 
in the space of snapshots: 

^ Af , (11) 

where 

= [«mn] = / = Cap^^snap 

J uj 

and 

5'°® = [S^„] = [ = i?^ap'5-Rsnap- 

J UJ 

The above integrals take into account the fracture distributions (see Eqn. (§)• We present the details of n 
later. 

Here A and S denote analogous fine-scale matrices as defined by 

Aij = / i<i[x)V(j)i • V(j)j dx Sij = / li{x)(j)i(j)j dx, 

J D J D 

where <f>i is the fine-scale basis function, and /5(x) is defined in the next subsection. 

To generate the offline space we then choose the smallest eigenvalues from Eqn. 0 and form the 
corresponding eigenvectors in the space of snapshots by setting (for /c = 1,..., M^), 

where are the coordinates of the vector 
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4.2 Global coupling 

In this subsection, we discuss the offline space and the variational formulation for a continuous Galerkin 
approximation of Eqn. 0 . We begin with an initial coarse space = span{x^}^]^. Denote N to be the 
number of coarse neighborhoods. Here, Xi standard multiscale partition of unity functions defined 

by 


— div (/i:(x) Vxi) = 0 K e uJi (12) 

Xi = gi on dK, 

for all K e uoi^ where gi is a continuous function on dK and is linear on each edge of dK. We note that 
pointwise energy H required for the eigenvalue problems is defined as 

N 

i=l 


where H denotes the coarse mesh size. 

We then multiply the partition of unity functions by the eigenfunctions in the offline space to 
construct the resulting basis functions 

^i,k = 1 <i < N and 1 < k < (13) 

where denotes the number of offline eigenvectors that are chosen for each coarse node i. We note that 
the construction in Eqn. ([T^ yields continuous basis functions due to the multiplication of offline eigenvectors 
with the initial (continuous) partition of unity function. Next, we define the continuous Galerkin spectral 
multiscale space as 

Wff = span{'0i^/c : 1 <i < N and 1 < A: < M^}. (14) 

Using a single index notation, we may write Uoff = where Nc = Ylf=i denotes the total 

number of basis functions in the space Uoff- We also construct an operator matrix Rq = ..., (where 

ijji are used to denote the nodal values of each basis function defined on the fine grid), for later use in this 
subsection. 

Below, we will display in detail the multiscale formulation using DEM and EEM, respectively. 


4.2.1 Multiscale DFM approach 

We seek i^ms(^) = G Kff such that 

a{ums^v) = {f,v) for all v G Uoff. (15) 


Note that the offline space 14ff is an approximation of all the nodal basis, 
fractures. Therefore, the variational form in (15) yields the following linear 


including the ones on the long 
algebraic system 


= Fo, 


(16) 


where Uq denotes the nodal values of the discrete solution, and 

Ao = [aij] = / k{x) dx and Fq = [//] = / fi^i dx. 

J D J D 


Using the operator matrix F^, we may write = RqARq and Fq = FqF, where A and F are the standard, 
hue-scale stiffness matrix and forcing vector corresponding to the form in Eqn. 0. We also note that the 
operator matrix may be analogously used in order to project coarse scale solutions onto the hne grid. 
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4.2.2 Multiscale EFM approach 

Similar to the multiscale DFM approach proposed above, we seek Ums(^) = G Kff such that 

a{ums,v) = {f,v) for all v G Vff- 

However, the nodal basis on the long fractures are excluded from the offline space Vqs, i.e., the fine-scale 
nodal basis are used on the long fractures. It follows from Eqn. 


RoAjnRo RoBlr 
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jNf _ 


(17) 


Note that the operator matrix may be analogously used in order to project coarse scale solutions onto the 
fine grid. 


5 Numerical result 

In this section, we will present several numerical experiments to show the performance of multiscale DFM 
and multiscale EFM approaches in the fracture modeling. The simulation results using multiscale DFM is 
shown in Subsection |5.1[ In Subsection |5.2[ a few tests are conducted to verify the performance of multiscale 
EFM approach. Moreover, we combine those two approaches in the modeling of short and long fractures 
and list the results in Subsection |5.3l Further, the relation between the basis selection and the types of 

We also present a randomized snapshot calculations to reduce the 
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fractures is studied in Subsection 
computational cost associated with computing the snapshot space. 

We take the domain D to be a square, set the forcing term / = 0 and impose a bilinear boundary condition 
for the problem Q. In our numerical simulations, we use a 10 x 10 coarse grid, and each coarse grid block 
is further divided into 10 x 10 fine grid blocks. Thus, the whole computational domain is partitioned into 
a 100 X 100 fine grid. The hue-scale solution is obtained by discretizing 0 following DFM with piecewise 
bilinear elements on the hne grid and linear one-dimensional elements over the fractures. 

We recall that Wff denotes the offline space; i^snap and denote the hne-scale, snapshot and offline 
solutions, respectively. In the tables below, we will compute the error u — using the relative error 
and the energy relative error, which are dehned as 


\u - UoE\\lI{D) • = 


11'^ — '^off||L2(y) 

ll'^IU2(y) 




Hl{D) 


i{u - Uofi.u- r^off )2 
a{u, u)^ 


(18) 


where the weighted L^-norm is dehned as ||i^||L 2 (y) = We will also compute the error i^snap~'^off 

using the same norms 


'‘'snap 


■ '^off||L2(D) 


ll'^snap - Moff||L=^(V) 
ll'^snap \\l^(V) 


^snap '^off 


a(u 


snap 




<3'(rf-snap5 '^snap) ^ 


(19) 


5.1 Numerical results with DFM 


In Fig. we depict three fracture helds used in the simulations below. In Fi g. |4(a)[ there is at most one 
fracture in each coarse block without crossing the coarse edges, while in Figs. |4(b) and 4(c)[ the fractures 
are more complicated and intersect the edges of coarse blocks. The simulation result with fractures in 
Fig. g ^with 1 basis per coarse node is very good (with the energy error of 2.52%) because their ehects 
can be localized. The simulation results with the other two fracture systems are shown in Tables and 
respectively. We notice that one basis function per node does not give a satisfactory for more complex 
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fracture systems. We also show the fine-scale solutions and the snapshot solutions (i.e. the coarse-scale 
solution corresponding to the largest offline space we select) corresponding to fractures in Fig. which are 
shown in Figs. [^and|^ respectively. 

The convergence history of fracture fields in Figs. |4(b)| and 4(c)| are displayed in Tables and We 
observe that the offline solution will converge to the fine-scale solution as we enrich the offline space in each 
coarse neighborhood. The energy error decreases from 27.25% to 7.4% as we add 4 more local basis in each 
coarse neighborhood from Table In Tables and we show the numerical results for fracture fields in 
Figs. |5(a)| and |5(b)[ 
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(a) Fracture 1 


(b) Fracture 2 


(c) Fracture 3 


Figure 4: Permeability fields with fractures. Fractures are shown with red, while the background is blue. 
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(a) Fracture 1 


(b) Fracture 2 


(c) Fracture 3 


Figure 5: Permeability fields with fractures. Fractures are shown with red, while the matrix is the blue 
region. 
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(a) Fine-scale solution 1 



(b) Fine-scale solution 2 
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(c) Fine-scale solution 3 


Figure 6: Fine-scale solutions corresponding to the permeability fields in 


Fig-i respectively. 
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(a) Fine-scale solution 1 (b) Fine-scale solution 2 (c) pressure solution 

Figure 7: Fine-scale solutions corresponding to the permeability fields in Fig. respectively. 



(a) Coarse-scale solution 1 (b) Coarse-scale solution 2 (c) Coarse-scale solution 3 

Figure 8: Coarse-scale solutions corresponding to the permeability fields in Fig. respectively 


We have also performed a few simulations with uniform triangular meshes with the fracture field depicted 
in Fig. |5(c)| We depict the fine-scale velocity components in Fig. The coarse-scale velocity component 
of this fracture field is similar to those fracture fields above, and thus are not presented. Triangular meshes 
can also be used in the simulation. 


5.2 Numerical results with EFM 


In this subsection, the numerical results using EFM are presented. As mentioned above, the main advantage 
of this model is that the fracture does not need to align with fine grid boundaries. The numerical results 
are shown in Figs. 10 and 11 Note that for this model, we have unique nodes (independent of the matrix 


nodes) for the fractures. The number of long fractures are quite few (e.g. there is only one fracture in Fig. 
Tol). Therefore, it is not necessary for the upscaling over the long fractures. The classical multiscale finite 


dim(Kff) 

h - «off|| (%) 

ll'^snap '^offll (%) 

L'UD) 

H^{D) 

Ll{D) 

mm 

121 

2.45 

27.25 

2.38 

26.13 

202 

0.63 

14.43 

0.54 

12.29 

283 

0.36 

11.74 

0.28 

9.01 

364 

0.29 

10.61 

0.21 

7.49 

445 

0.16 

7.40 

— 

— 


Table 1: Convergence history for GMsFEM with different coarse spaces dimensions corresponding to the 
permeability field in Eig. |4(b)[ 
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dim(yoff) 

h - «off|| (%) 

ll'^snap '^offll (%) 

LUD) 


L'UD) 

miD) 

121 

0.58 

11.77 

0.46 

9.56 

202 

0.28 

8.93 

0.17 

5.74 

283 

0.24 

7.98 

0.10 

4.11 

364 

0.17 

6.84 

— 

— 


Table 2: Convergence history for GMsFEM with different coarse spaces dimensions corresponding to the 
permeability field in Fig. |4(c)[ 


dim(14ff) 

h - «off|| (%) 

ll'^snap '^offll (^) 

L'UD) 

H^{D) 

Ll{D) 

mm 

121 

1.61 

24.46 

1.53 

23.21 

202 

0.45 

13.33 

0.37 

10.98 

283 

0.27 

10.33 

0.17 

7.07 

364 

0.20 

8.58 

0.09 

4.13 

445 

0.17 

7.51 

— 

— 


Table 3: Convergence history for GMsFEM with different coarse spaces dimensions corresponding to the 
permeability field in Fig. |5(a)| 


element method is applied in Fig. 


error and errors are 0.3327% and 0.1734%. 
field. Similarly, for the simulation depicted in Fig. 
and 0.1409% for energy error and error respectively. 


10 to obtain the coarse-scale solution. The corresponding weighted 

Therefore, no spectral problem is needed for this fracture 
the GMsFEM algorithm gives an error of 0.1344% 
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5.3 Numerical results for mutiscale DFM and EFM approach 


In this subsection, we will display several numerical results using DFM and EFM jointly. A heterogeneous 
medium with different lengths of fractures is illustrated in Fig 12(a)| The coarse-scale solution and the fine- 
scale solution are shown in Figs |12(b)| and 12(c)[ respectively. The fine-scale solution is calculated by the 
DFM over the fine mesh. We observe that the coarse-scale solution is a good approximation of the fine-scale 
solution. Besides, the relative energy error and relative errors are shown in Table We see that the 
energy relative error and relative error are 9.17% and 0.02% when the dimension of the coarse system is 


364. 

We test another heterogeneous fractured medium in Fig. |13(a)| Gompared with the previous heteroge¬ 
neous medium, this one has a longer fracture with a curved shape. The coarse-scale solution and fine-scale 
solution are shown in Figs. |13(b)| and |13(c)| We present the errors in Table 


diin(14ff) 

(%) 

ll'^snap '^offll (%) 

Lum 

mm 

mm 

mm 

121 

1.28 

19.54 

1.22 

18.57 

202 

0.25 

9.39 

0.20 

7.23 

283 

0.19 

8.00 

0.13 

5.31 

364 

0.16 

7.10 

0.09 

3.82 

445 

0.12 

5.98 

— 

— 


Table 4: Gonvergence history for GMsFEM with different coarse spaces dimensions corresponding to the 
permeability field in Fig. |5(b)[ 
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Figure 9: Velocity field corresponding to the permeability field in Fig. |5(c)[ 
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(a) Fracture field 


(b) Fine-scale solution 


(c) coarse-scale solution 


Figure 10: Numerical results using GMsFEM with EFM. 

5.4 Adaptive method 

In this subsection, we will discuss adaptive strategies for generating multiscale basis functions efficiently. 
Eirst, we will not use extra basis functions in the regions with long fractures. This intuitive approach is 
based on previous findings uni. Next, we will use the error indicators that we proposed in [4], which identify 
the regions adaptive enrichment is needed. 

We consider the fracture fields shown in Eigs. |4(b)| and 4(c) As shown in these figures, both fields have 



(a) Fracture field 
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(b) Fine-scale solution (c) coarse-scale solution 


Eigure 11: Numerical results using GMsEEM with EME. 
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Figure 12: Permeability field (left), coarse-scale solution (middle) and fine-scale solution (right) using the 
DFM model for short fractures and EFM for long fractures. 


dim(Uoff) 

h-«off|| (%) 

ll'^snap '^offll (^) 

LUD) 


L'UD) 

miD) 

121 

3.14 

96.16 

3.05 

85.15 

202 

2.22 

73.10 

2.15 

63.04 

283 

0.04 

12.59 

0.03 

5.66 

364 

0.02 

9.17 

0.01 

2.41 

445 

0.01 

6.66 

— 

— 


Table 5: Convergence history for GMsFEM with different coarse spaces dimensions corresponding to the 
permeability field in Eig. using DEM for resolving the short fractures and hierarchical method for the 
long fractures. 


channels and isolated inclusions in certain coarse neighborhoods. The local coarse spaces contain enough 
information about isolated inclusions because of the usage of multiscale partition of unity functions in the 
construction of global generalized multiscale basis. Therefore, the adaptive enrichment process takes place 
in the regions with long channels. 

We begin with 1 basis in each coarse node and examine the fracture field shown in Eig |4(b)| following 
the reasoning above. We conclude that the coarse nodes with adaptive enrichment are (i, j), 3 < i < 9 and 
5 < J < 9. Adding 3 more basis to those coarse nodes, we can get the energy error of 8.55% with the coarse 
space dimension of 226. Compare with results listed in Table (coarse space of dimension 364 and the 
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(a) Fracture system 


(b) Coarse-scale solution 


(c) Fine-scale solution 


Eigure 13: Permeability field (left), coarse-scale solution (middle) and fine-scale solution (right) using the 
DEM model for short fractures and EEM for long fractures. 
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dim(yoff) 

h - «off|| (%) 

ll'^snap '^offll (%) 

LUD) 


L'UD) 

miD) 

121 

0.63 

37.85 

0.54 

28.69 

202 

0.52 

32.23 

0.45 

22.55 

283 

0.03 

10.15 

0.01 

2.54 

364 

0.02 

8.13 

0.0008 

0.09 

445 

0.02 

7.55 

— 

— 


Table 6: Convergence history for GMsFEM with different coarse spaces dimensions corresponding to the 
permeability field in Fig. [fusing DFM for resolving the short fractures and EFM for long fractures. 


energy error of 10.61%), we observe that with this adaptive enrichment, we can get a better offline solution 
with a smaller solution space. 

We also test the field in Fig. |4(c)[ The coarse nodes with long fractures are (i, j), 3 < i < 9 and 
6 < i < 8. Adding 2 more basis to those coarse nodes, we can get the energy error of 7.70% with the coarse 
space dimension of 163. Compare with results listed in Tablethe energy error is 7.98% with a coarse space 
of dimension 283. We can see that with this adaptive algorithm, we can get a better offline solution with a 
smaller solution space. 

Next, we use the i7“^-based error indicator proposed in |4]. We start with 1 basis per coarse node and 
take 0 = 0.7 in the adaptive algorithm proposed in [4]. Comparing Table with Tabl e we notice that 
there is a slight gain regarding to the energy error for this adaptive algorithm. In Table [SJthe energy error 
is 4.51% for an offline space of dimension 445, while 4.39% for 424 in the adaptive algorithm. 

Next, we study the adaptive algorithm for the field in Fig. |4(c)| From results in Table and Table 
we conclude that the adaptive algorithm can improve the performance of GMsFEM. Using an offline space 
with 80 less basis, the adaptive enrichment algorithm can still have a better performance. 


diin(Vff) 

llti-Refill (%) 

ll'^snap '^offll (%) 

LUD) 


L'iiD) 


121 

1.61 

24.46 

1.53 

23.21 

159 

0.44 

13.49 

0.36 

11.15 

230 

0.27 

10.36 

0.20 

7.32 

362 

0.18 

8.37 

0.14 

5.50 

424 

0.15 

7.45 

— 

— 


Table 7: Convergence history for GMsFEM with different coarse spaces dimensions using adaptive algorithm 
corresponding to permeability field in Fig. |5(a)[ 


dim(Kff) 

h - «off|| (%) 

ll'^snap '^offll (%) 

L'iiD) 

H^{D) 

L'iiD) 


121 

1.28 

19.54 

1.22 

18.57 

169 

0.25 

9.35 

0.20 

7.14 

259 

0.19 

8.01 

0.13 

5.27 

310 

0.14 

6.85 

0.09 

3.72 

370 

0.10 

5.49 

— 

— 


Table 8: Gonvergence history for GMsFEM with different coarse spaces dimensions using adaptive algorithm 
corresponding to permeability field in Fig. |5(b)[ 
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Figure 14: Illustration of a coarse neighborhood and oversampled domain. Here, K is a coarse-grid block, 
uji is a coarse neighborhood of and uj^' is an oversampled region 

5.5 Randomized snapshots 

In this section, we will investigate the oversampling randomized algorithm proposed in [2]. The main 
advantage of this algorithm is that it uses much fewer snapshot basis functions to calculate the offline space. 
Besides, the oversampling strategy is used to reduce the sub grid errors due to boundary conditions imposed 
to obtain the snapshot basis. We refer to Fig. [^for an illustration of oversampling domain. 

In our simulation, we set the oversampling size t = 2 (two fine-grid layer around the coarse region) and 
buffer number = 2 for each coarse neighborhood uJi. I.e., we generate n + 2 snapshot functions in each 
coarse neighborhood when computing n basis functions. For the sake of completeness, we list the algorithm 
here in Table [9l 


Table 9: Randomized GMsFEM Algorithm ([2]) 


Input: 

output: 

1 . 

2 . 

3. 

4. 


Fine grid size h, coarse grid size oversampling size t, buffer number for each uji^ 
the number of local basis functions for each 
Coarse-scale solution uh- 

Generate oversampling region for each coarse block: T^, T^, and ujf ; 

Generate random vectors ri and obtain randomized snapshots in 00 ^'; 

Add a snapshot that represents the constant function on ; 

Obtain k^^ offline basis by a spectral decomposition restricted to the original domain 
Gonstruct multiscale basis functions and solve it. 


The numerical results using the permeability field shown in Fig. |I3(a)| are presented in Table 10 In this 
table, the first column shows the dimension of the offline space, and the second columns displays the ratio 
between the number of the randomized snapshot basis and the full snapshot basis. Then, in the following two 
columns, the weighted and energy error using the full snapshots are shown. At last, the results using the 
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randomized snapshots alone is listed in the last two column. Comparing these results with those using the 
full snapshot space, we deduce that the multiscale randomized algorithm performs well. The relative energy 
error is 15.82% using the randomized snapshots, which account to 6.25% of the full snapshots. However, we 
can only get 13.50% relative energy error if the full snapshots are used instead. 


Table 10: Numerical results using randomized snapshots, pbf = 4. 


dim(yoff) 

Snapshot ratio (%) 

All snapshots (%) 

Randomized snapshots (%) 

LiiD) 

Hk{D) 

LUD) 

h:{d) 

283 

5.21 

7.31 

15.55 

8.07 

16.85 

364 

6.25 

5.05 

13.50 

6.60 

15.82 

445 

7.29 

4.46 

12.59 

5.80 

14.48 


6 Discussions and Conclusions 

In this paper, we develop a multiscale finite element method for flows in fractured media. Our approach 
is based on Generalized Multiscale Finite Element Method (GMsFEM). Multiscale basis functions are con¬ 
structed in the offline stage via local spectral problems. To represent the fractures on the fine grid, we 
consider two approaches (1) Discrete Fracture Model (DEM) (2) Embedded Fracture Model (EFM) and 
their combination. The proposed procedure automatically selects multiscale basis functions via local spec¬ 
tral problems. 

Numerical results are presented where we demonstrate how one can adaptively add basis functions in 
the regions of interest based on error indicators. We consider various cases with long and short fractures. 
Our numerical results show that GMsFEM with much fewer degrees of freedom compared to the fine-scale 
simulations can capture the solution behavior accurately. Furturemore. we discuss the use of randomized 
snapshots ([2]) which substantially reduces the offline computational cost. 
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